function [values2,ofv,tsub,infeas,binom,p,nn,timertime] = UP6c_p()

%load remaining data
data_struct = load('data6c.mat');
segNP = double(data_struct.segNP);
Y = data_struct.Y;
cap = double(data_struct.cap);
eps = data_struct.eps;
bound = data_struct.bound;
upMat = data_struct.upMat;
os = double(data_struct.os);
uniqSeg = double(data_struct.uniqSeg);
nyp = double(data_struct.nyp);
prices = double(data_struct.prices);
capcon = double(data_struct.capcon);
segScale = double(data_struct.segScale);
segQ = double(data_struct.segQ);
rem_cap = double(data_struct.rem_cap);
uniqSeg_size = size(uniqSeg);
nNewSegments = uniqSeg_size(1);
prices_size = size(prices);
nProducts = prices_size(1);
[y1,y2,y3]=ind2sub(size(upMat),find(upMat));
nRes = length(rem_cap);
inst = double(data_struct.inst);

prices_lb = zeros(nProducts,nProducts,nNewSegments);
for i=1:nNewSegments
     prices_lb(y1(i),y2(i),i) = prices(y1(i));
end

u = zeros(nNewSegments,3+nProducts);
usub = zeros(nNewSegments,nProducts);
for i = 1:nNewSegments
    u(i,:) = [uniqSeg(i,1)+1,y1(i),y2(i),os(uniqSeg(i,3)+1,:)];
    usub(i,:) = os(uniqSeg(i,3)+1,:);
    usub(i,y2(i)) = 0;
end

% Matrix A
anzNBcap = find(rem_cap~=0,1,'last')-1;
anzNByp = nNewSegments;
anzEV = 2*nNewSegments;
A = zeros(anzNBcap+anzNByp,anzEV);
for i = 1:anzNBcap
    h = i+1;
    for n = 1:nNewSegments
        A(i,n) = (capcon(y1(n),i)-capcon(y1(n),h));
    end
end
A(anzNBcap+1:anzNBcap+anzNByp,1:nNewSegments) = diag(ones(nNewSegments,1));
A(anzNBcap+1:anzNBcap+anzNByp,nNewSegments+1:2*nNewSegments) = -diag(Y(sub2ind(size(Y),y1,y2,y3)));

%b, Aeq, beq
b = [rem_cap(2:anzNBcap+1)';zeros(anzNByp,1)];
Aeq = [];
beq = [];

%OF
of = @of6c_p;

[~,q2] = find(capcon~=0);
if bound == 'lb'
    x0 = [zeros(nNewSegments,1);0.1*ones(nNewSegments,1)];
elseif bound == 'ub'
    x0 = [zeros(nNewSegments,1);(1-0.1)*ones(nNewSegments,1)]; 
    x0(nNewSegments) = min(Y(sub2ind(size(Y),y1(nNewSegments),y2(nNewSegments),nNewSegments)),rem_cap(q2(y2(nNewSegments))));
elseif bound == 'mb'
    x0 = [zeros(nNewSegments,1);(1-0.1)*ones(nNewSegments,1)];
    first = find(y2==max(y2),1,'first');
    x0(first) = min(Y(sub2ind(size(Y),y1(first),y2(first),first)),rem_cap(q2(y2(first))));
end

lb = [zeros(nNewSegments,1);eps*ones(nNewSegments,1)];
ub = [Y(sub2ind(size(Y),y1,y2,y3));(1-eps)*ones(nNewSegments,1)];

% fmincon
options = optimoptions('fmincon','MaxIterations',5000,'MaxFunctionEvaluations',30000);
tic;
[x,ofv,ef] = fmincon(of,x0,A,b,Aeq,beq,lb,ub,[],options);
tsub = toc;

nn = x(1:nNewSegments);
p = x(nNewSegments+1:end);
values2 = zeros(nNewSegments,6);
values2(:,5) = x(nNewSegments+1:end);
values2(:,6) = Y(sub2ind(size(Y),y1,y2,y3));
for i=1:nNewSegments
    values2(i,1:2) = [i,y2(i)];
    values2(i,4) = x(i);
    if nProducts == 6
        values2(i,3) = segScale(u(i,1))*(reallog(1-x(nNewSegments+i))...
                        -reallog((exp(segNP(u(i,1)))+u(i,4)*exp(1/segScale(u(i,1))*(segQ(1,u(i,1))-prices(1)))...
                               +u(i,5)*exp(1/segScale(u(i,1))*(segQ(2,u(i,1))-prices(2)))...
                               +u(i,6)*exp(1/segScale(u(i,1))*(segQ(3,u(i,1))-prices(3)))...
                               +u(i,7)*exp(1/segScale(u(i,1))*(segQ(4,u(i,1))-prices(4)))...
                               +u(i,8)*exp(1/segScale(u(i,1))*(segQ(5,u(i,1))-prices(5)))...
                               +u(i,9)*exp(1/segScale(u(i,1))*(segQ(6,u(i,1))-prices(6))))...
                        -(1-x(nNewSegments+i))*(exp(segNP(u(i,1)))+usub(i,1)*exp(1/segScale(u(i,1))*(segQ(1,u(i,1))-prices(1)))...
                               +usub(i,2)*exp(1/segScale(u(i,1))*(segQ(2,u(i,1))-prices(2)))...
                               +usub(i,3)*exp(1/segScale(u(i,1))*(segQ(3,u(i,1))-prices(3)))...
                               +usub(i,4)*exp(1/segScale(u(i,1))*(segQ(4,u(i,1))-prices(4)))...
                               +usub(i,5)*exp(1/segScale(u(i,1))*(segQ(5,u(i,1))-prices(5)))...
                               +usub(i,6)*exp(1/segScale(u(i,1))*(segQ(6,u(i,1))-prices(6))))))...
                               +segQ(u(i,3),u(i,1));
   elseif nProducts == 12
        values2(i,3) = prices_lb(y1(i),y2(i),i)...
                        +segScale(u(i,1))*(reallog(1-x(nNewSegments+i))...
                        -reallog(exp(segNP(u(i,1)))+u(i,4)*exp(1/segScale(u(i,1))*(segQ(1,u(i,1))-prices(1)))...
                               +u(i,5)*exp(1/segScale(u(i,1))*(segQ(2,u(i,1))-prices(2)))...
                               +u(i,6)*exp(1/segScale(u(i,1))*(segQ(3,u(i,1))-prices(3)))...
                               +u(i,7)*exp(1/segScale(u(i,1))*(segQ(4,u(i,1))-prices(4)))...
                               +u(i,8)*exp(1/segScale(u(i,1))*(segQ(5,u(i,1))-prices(5)))...
                               +u(i,9)*exp(1/segScale(u(i,1))*(segQ(6,u(i,1))-prices(6)))...
                               +u(i,10)*exp(1/segScale(u(i,1))*(segQ(7,u(i,1))-prices(7)))...
                               +u(i,11)*exp(1/segScale(u(i,1))*(segQ(8,u(i,1))-prices(8)))...
                               +u(i,12)*exp(1/segScale(u(i,1))*(segQ(9,u(i,1))-prices(9)))...
                               +u(i,13)*exp(1/segScale(u(i,1))*(segQ(10,u(i,1))-prices(10)))...
                               +u(i,14)*exp(1/segScale(u(i,1))*(segQ(11,u(i,1))-prices(11)))...
                               +u(i,15)*exp(1/segScale(u(i,1))*(segQ(12,u(i,1))-prices(12)))...
                        -(1-x(nNewSegments+i))*(exp(segNP(u(i,1)))+usub(i,1)*exp(1/segScale(u(i,1))*(segQ(1,u(i,1))-prices(1)))...
                               +usub(i,2)*exp(1/segScale(u(i,1))*(segQ(2,u(i,1))-prices(2)))...
                               +usub(i,3)*exp(1/segScale(u(i,1))*(segQ(3,u(i,1))-prices(3)))...
                               +usub(i,4)*exp(1/segScale(u(i,1))*(segQ(4,u(i,1))-prices(4)))...
                               +usub(i,5)*exp(1/segScale(u(i,1))*(segQ(5,u(i,1))-prices(5)))...
                               +usub(i,6)*exp(1/segScale(u(i,1))*(segQ(6,u(i,1))-prices(6)))...
                               +usub(i,7)*exp(1/segScale(u(i,1))*(segQ(7,u(i,1))-prices(7)))...
                               +usub(i,8)*exp(1/segScale(u(i,1))*(segQ(8,u(i,1))-prices(8)))...
                               +usub(i,9)*exp(1/segScale(u(i,1))*(segQ(9,u(i,1))-prices(9)))...
                               +usub(i,10)*exp(1/segScale(u(i,1))*(segQ(10,u(i,1))-prices(10)))...
                               +usub(i,11)*exp(1/segScale(u(i,1))*(segQ(11,u(i,1))-prices(11)))...
                               +usub(i,12)*exp(1/segScale(u(i,1))*(segQ(12,u(i,1))-prices(12))))))...
                               +segQ(u(i,3),u(i,1))-prices(u(i,2));
    elseif nProducts == 2
        values2(i,3) = prices_lb(y1(i),y2(i),i)...
                    +segScale(u(i,1))*(log(1-x(nNewSegments+i))...
                    -reallog(exp(segNP(u(i,1)))+u(i,4)*exp(1/segScale(u(i,1))*(segQ(1,u(i,1))-prices(1)))...
                           +u(i,5)*exp(1/segScale(u(i,1))*(segQ(2,u(i,1))-prices(2)))...
                    -(1-x(nNewSegments+i))*(exp(segNP(u(i,1)))+usub(i,1)*exp(1/segScale(u(i,1))*(segQ(1,u(i,1))-prices(1)))...
                           +usub(i,2)*exp(1/segScale(u(i,1))*(segQ(2,u(i,1))-prices(2))))))...
                           +segQ(u(i,3),u(i,1))-prices(u(i,2));
    end
end

infeas = ef;

values2 = values2(:,1:4); 
binom = 0;
 